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-of  temperature  can  be  used  to  calculate  whether  or  not  a  given  energy  source  is  adequate  to  ignite 
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mixture  of  H2-02-N2"for  two  values  of  R0  which  show  that  a  model  must  be  constructed  for  a 
quench  volume  in  order  to  complete  the  similarity  solution  calibration. 
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THEORETICAL  AND  COMPUTATIONAL  APPROACH  TO 
MODELLING  FLAME  IGNITION 

Introduction 

Ignition  of  fuel-oxidizer  mixtures  occur  when  an  external  source 
of  energy  initiates  interactions  among  the  controlling  convective, 
transport  and  chemical  processes.  Whether  the  process  results  in 
deflagration,  detonation,  or  is  simply  quenched  depends  on  the  in¬ 
tensity,  duration,  and  volume  affected  ny  an  external  heat  source. 

Ignition  also  will  depend  on  the  initial  ambient  properties  of  the 
mixture  which  determine  the  chemical  induction  time  and  the  heat  release 
per  gram  of  material.  This  ignition  is  a  complicated  phenomena 
whose  occurrence  for  a  specific  mixture  of  fuel  and  oxidizer  depends 
strongly  on  diffusive  and  chemical  parameters  which  are  often  very 
poorly  known.  A  convenient,  inexpensive  way  to  estimate  whether  a 
mixture  will  ignite  given  a  heat  source  intensity,  duration,  and  volume 
would  be  both  a  valuable  laboratory  tool  and  a  useful  learning  device. 

This  paper  represents  a  summary  of  the  work  to  date  on  the 
theoretical  and  computational  effort  at  the  Naval  Research  Laboratory 
in  the  study  of  the  ignition  of  homogeneous  gas  phase  mixtures.  First 
the  essential  properties  of  a  detailed,  time-dependent,  numerical  model 
are  presented.  Then  a  description  is  given  of  a  closed  form 
similarity  solution  for  the  nonlinear  time-dependent  slow  flow  equations 
which  forms  the  basis  for  a  simple,  time-dependent  "hand-calculator" 
model  of  localized  ignition.  The  similarity  model  requires  minimal 
chemical  and  physical  input  and  contains  two  constants  which  must 
be  calibrated:  the  radii,  or  fraction  of  the  time-dependent  similarity 
solution  radius, at  which  the  thermal  conductivity  and  induction  paramenters 
ere  evaluated.  Finally,  camparisons  between  two  models  are  described 

for  ignition  of  an  H2_02_N2  mixture  *  The  composition  of  the 
Note:  Manuscript  submitted  October  9,  1979. 


1 


material  and  the  length  of  time  of  energy  deposition  are  held  constant 


but  the  total  energy  deposition  and  the  radius  of  deposition  are  allowed 
to  vary. 

The  Detailed  Simulation  Model 

In  this  section  we  describe  the  detailed  numerical  flame  model  and 
present  a  sample  of  calculations  performed  with  it.  The  model  permits 
a  wide  variety  of  geometric , initial ,  boundary  and  time-varying  energy 
input  conditions  and  was  specifically  developed  to  study  the  various 
physical  and  chemical  processes  which  control  flame  initiation  and 
quenching.  Whereas  propagation  of  a  laminar  flame  has  been  successfully 
treated  by  steady  state  models,  the  details  of  initiation  and  quenching 
are  inherently  time-dependent  phenomena  and  require  detailed  time-dependent 
models. 

The  model  we  have  developed  solves  in  one  dimension  the  full  set 
coupled  conservation  equations  for  total  mass  p,  monentum  pv,  and  energy 
E  as  well  as  the  individual  species  densities  {n_,}  : 

~  =  -  V  •  pv  (1) 

Ot 


V  •  n.  V.  -V  •  n,v  +  P.  -  Q.nj 
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The  quantity  v  is  the  fluid  velocity,  the  {v^.}  are  the  diffusion  velocities, 

and  {P  }  and  {Q.}  refer  to  the  chemical  production  and  loss  processes 
j  1 

for  the  individual  species  j .  The  quantities  n  and  1  are  the  mixture 

m  m 

viscosities  or  thermal  conductivities  using  the  expressions  of  Wilke 

(1950)  and  Mason  and  Saxena,  (1958) .  The  superscript  "T"  in  the  last  term 

of  Eq.  (3)  indicates  that  the  transpose  is  to  be  taken.  The  quantity 

Qq  represents  the  contributions  to  the  thermal  conductivity  due  to 

diffusion  processes  and  it  is  a  function  of  {V^ } ,  the  binary  molecular 

diffusion  coefficients  {D  },  the  thermal  diffusion  coefficients  {D}}, 

jk 

and  the  enthalpies  {h^}.  Solution  of  the  diffusion  velocities  includes 

the  binary  and  thermal  diffusion  terms  as  well  as  the  pressure  gradient 

term, which  in  this  flame  is  negligible.  We  also  assume  here  that  we  are 

dealing  with  mixtures  of  ideal  gases  so  that  the  pressure,  P,  may  be  written 

P  =  Nk  T  where  N  is  the  total  number  density,  k  is  Boltzmann's  constant 
3  B 

and  T  is  the  temperature.  The  model,  however,  is  not  restricted  to  ideal 
gases  and  in  fact  any  equation  of  state  may  be  used.  Finally,  we  assume 
that  these  are  no  radiative  losses  in  the  system. 

The  technique  for  solving  Eqs  (1)-  (4)  is  based  on  the  method 
of  asymptotic  timestep  splitting.  The  individual  chemical  and  physical 
processes  are  integrated  separately  by  the  fastest  and  most  accurate 
algorithms  and  then  asympotically  coupled  together.  Advantages  of 
this  procedure  over  the  more  standard  global  implicit  approaches 
(e.g.  Westbrook,  1978;  Lund,  1978)  include  faster  calculations  for 
modest  and  large  numbers  of  species  and  more  flexibility  in  problem 
configuration.  The  work  presented  here  shows  quite  convincingly  that 
our  asymptotic  coupling  mehtods  work  accurately,  reliably,  and 
efficiently  despite  mistaken  grumbling  to  the  contrary. 
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The  convective  transport  is  solved  by  the  algorithm  ADINC  (  Boris, 

1979), an  implicit  Lagrangian  hydrodynamics  algorithm  for  subsonic  flows. 

The  method  allows  for  an  arbitrary  equation  of  state  and  gives  an  accurate 
representation  of  material  interfaces.  Thus  it  is  also  well  suited  to 
the  analysis  of  heterogeneous  combustion.  ADINC  communicates  compression  and 
expansion  accross  the  system  implicity  and  yet  maintains  the  steep 
qradients  in  species  and  temperature  required  for  combustion  modelling. 

The  implicit  pressure  calculation  overcomes  the  Courant  timestep 
limit.  Tests  of  the  full  detailed  flame  model  in  the  limit  of  no 
chemical  reactions  of  diffusive  transport  are  essentially  tests  of  the  ADINC 
algorithm.  A  number  of  such  tests  have  been  documented  in  Boris,  1979. 

We  further  note  that  this  algorithum  lends  itself  to  an  adaptive 
gridding  method  in  which  cells  are  inserted  or  deleted  according 
to  externally  specified  physical  conditions  in  the  flow. 

The  chemical  interactions  are  described  by  a  set  of  nonlinear 
coupled  ordinary  differential  equations.  For  this  ignition  calibration, 
we  have  used  the  H2-02  reacti°n  scheme  given  in  Table  I  which  involves 
the  eight  reactive  species,  0^,  0,  H,  OH,  H02,  ^2°2'  H2°'  and  diluent 
which  is  chosen  to  be  •  The  thermochemical  properties  of  these  species 
were  taken  from  the  JANAF  tables  (Stull  and  Prophet,  1971) . 

The  final  chemical  rate  scheme  was  arrived  at  after  an  extensive  study 
of  the  literature  especially  aimed  at  updating  the  work  of  Baulch  (1972) . 

The  mechanism  was  then  tested  extensively  using  the  intergration  method 
VSAIM,  which  is  a  fully  vectorized  version  of  the  selected  asymtotic 
intergration  method,  CHEMEQ,  developed  by  Toung  (1977,1979) .  Many 
of  these  tests  were  performed  in  the  limit  of  no  convective  or  diffusive 
transport,  which  allowed  detailed  comparisons  with  measured  induction 
times  over  a  wide  range  of  temperatures,  pressures,  and  stoichiometries. 
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Typical  results  obtained  by  Burks  and  Oran  (1980)  are  shown  in  Fig  (1) , 


which  shows  that  we  have  obtained  good  agreement  without  any  reaction 
rate  adjustments.  Thus  we  have  confidence  in  both  the  reaction  scheme 
and  the  intergration  method. 

A  new  method  has  been  devised  to  solve  the  {V_. }  accurately  and 
without  matrix  inversion  (Boris  and  Oran,  1978) .  The  diffusion  equations 
may  be  written 


S  .= 
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M  n  ,n 
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Then  it  may  be  shown  that  V_.  may  be  written 
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j T,  the  diffusion  coefficient  of  species  j  through  the  background 
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provided  by  the  sum  of  all  the  other  species,  is  defined  by 
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and  the  matrix  elements  of  A  are  given  by 
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This  algorithm  is  of  0(M  )  and  is  vectorized.  Thus  it  is  substantially 
3 

faster  than  0  (II  )  matrix  inversions  when  more  than  four  of  five  species 
are  involved. 

The  detailed  model  described  above  is  currently  written  in  ANSI 
standard  Fortran  and  vectorized  for  operation  in  the  Texas  Instruments 
Advanced  Scientific  Computer.  Running  time  is  between  .05  and  .08 
seconds  for  a  timestep  which  includes  evaluations  of  the  full  set  of 
Eqs.  (1)  -  (11)  for  45  cells. 

Results  from  a  typical  flame  calculation  in  cartesian  geometry  are  shown 
in  Figs (2) ,  (3)  and  (4).  Figure  (2)  shows  a  typical  temperature 
profile  for  a  flame  initiated  by  a  quadratic  temperature  distribution 
at  the  onset  of  the  calculation.  Also  shown  on  the  same  figure  is 
a  graph  of  the  cell  size,  Ax,  which  shows  that  the  resolution 
decreases  as  the  absolute  magnitude  of  temperature  gradient  increases. 

Species  profiles  are  shown  in  Fig. (3)  and  details  of  the  flame  front  in 
Fig  (4) .  We  note  that  the  intermediate  species  and  the  temperature 


gradients  do  not  all  peak  at  the  same  location  and  the 

structure  of  the  flamefront  is  clearly  distingiushable.  This  will  be 

discussed  more  thoroughly  in  a  future  paper. 

The  Similarity  Solution  Model 

The  basic  similarity  solution  is  derived  from  the  slow  flow  (Jones 
and  Boris,  1977,  Boris  and  Oran,  1978)  approximate  and  is  predicated 
on  the  assumption  that  energy  addition  to  the  sy.<'.  am  is  slow  enough 
that  no  appreciable  energy  is  deposited  in  the  heated  region  by  shocks. 
Thus  the  system  is  characterized  by  (1)  flow  velocities  which  are  small 
compared  to  the  speed  of  sound,  and  (2)  an  essentially  constant  pressure 
field.  The  energy  and  velocity  equations  may  then  be  written  as 

2  2 

=  -  yP l_  •.  v  +  7  .  yNk  kAT  +  S(t)e  _k  (t)r*  (12) 

We  also  require  the  continuity  equation  which  is  written  in  the 
form 


i  =  _  v  • 

P  dp  - 


v  • 


(13) 


In  Eq.  (12) ,  Y  is  the  ratio  of  heat  capacities  C^/C^,  assumed  here 
to  be  constant,  and  <  is  a  function  of  the  mixture  thermal  conductivity. 


*  -  y-i  i 


(T)  . 


(14) 
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The  last  term  on  the  right  hand  side  of  Eq.  (12)  is  the  source  term. 

Proper  choice  of  S(t)  ensures  that  a  given  amount  of  energy,  E^,  is 

4it  3 

deposited  in  a  certian  volume,  —  R  ,  in  a  time,  t  .  It  is  the  choice 

3  o  o 

of  this  Gaussian  profile  which  allows  us  to  obtain  a  "closed"  form 

similarity  solution  given  below  in  spherical  coordinates. 

Using  the  assumption  that  dP/  as  0,  an  algebraic  equation 

dt 

is  written  for  V_  .  v  from  Eq.  (12)  which  when  combined  with  Eq.  (13) 
gives 


I®  =  £(t)  -  k2(t)r2  +  V.  K  7  T  (15, 

T  dt  P  ~  t  U 

(16) 

(17) 

where  T  and  p  are  the  background  temperature  and  pressure,  respectively. 
00  00 

Thus  the  nonlinear  slow-flow  equations  including  expansions  and 
contractions  of  the  flow  have  been  converted  into  a  single  equation 
which  is  linear  in  the  logarithm  of  the  temperature. 

The  total  energy  of  the  system  at  any  instant  is  the  sum  of  the 
internal  energy  and  the  work  performed  in  expanding  the  heated  region.  It 
may  be  written 

00 

E (t)  =  J  4nr2dr  £  i  -  o  (r,t) 

O  poo  (18) 


where  P  is  the  background  pressure.  The  solution  is  then 


T(r , t)  =  T  e 


A(t)e 
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2  2 
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V  P. 


(Y-l)  k  (t) 


4ttx 


- - -  F  (A  (t)  ) 

(Y-l)  kit) 


(19) 


which  defines  the  intergral  F(A(t)).  Differentiating  this  we  find  that 


dE  _  n3'2  S(t) 

dt  <Y“Dk^  (t)  (20) 

which  may  be  equated  to 

dE  =  3E  3k  +  3E  3A 

dt  ~  3k  3t  3a  3t  •  (21) 

Thus  a  consistency  condition  has  been  specified  on  the  rates  of 
change  of  the  amplitude,  A(t) ,  and  the  scale  size  k  1(t)  for  the  heated 
region. 

If  the  fluid  velocity  v  is  then  expanded  such  that 

v(r)  asv^  (t)  r  ,  (22) 

that  is,  only  the  linear  term  is  kept,  two  coupled  ordinary  differntial 
equations  for  k  and  A  may  be  obtained  from  Eqs.  (15)-  (21), 


dk 

dt 


-  kv. 


dA  =  S(t) 
dt  yP 


-  2Kk' 


-  6<k  A 


(23) 


(24) 
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and  we  also  find  that 


V  -  F'  (0)  -F’  (A)  2  AF'  (A)  -F (A) 

1  F (A)  +  ZKKB  F (A) 


(25) 


This  last  expression  for  v  ,  Eq.  (25)  is  basically  the  result  of 
envoking  energy  conservation  through  equating  Eqs.  (20)  and  (21)  and 
is  the  first  calibration  in  the  model. 

The  model  requires  one  further  definition  in  order  to  predict  ignition. 
A  curve  of  chemical  induction  time  as  a  function  of  temperature  must 
be  included  in  order  to  define  the  induction  parameter, 

Kt)  =  f  _ 

JQ  tc(T(r  ,t'))  •  f  (26) 

is 

Ignition  "occurs"  when  I(t)  *  1  in  this  model,  which  an  exact  result  in 
the  limit  of  large  heat  source  and  constant  temperature  near  the  center  of 
the  heated  region.  A  simple  analytic  expression  for  Tc (T)  depending 

on  three  constants  has  been  derived  and  can  be  calibrated  using  as 
few  as  three  distinct  values  of  at  different  temperatures.  These 
values  may  be  obtained  from  detailed  kinetic  calculations,  a  few 
measured  points  in  the  curve,  or  from  educated  guesses. 
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Model  Comparisons 

In  order  to  compare  the  results  from  the  detailed  similation  model 

to  that  of  the  similarity  model,  the  detailed  model  has  been  configured 

for  spherical  symmetry  with  open  boundary  conditions.  Energy  is  then 

deposited  in  a  Gaussian  profile  with  a  characteristic  radius  which  closely 

matches  the  one  derived  from  the  similarity  model.  Energy  deposition 

in  both  models  is  assumed  linear  in  time  at  a  rate  determined  by  requiring 

an  energy  Eq  to  be  deposited  in  a  time  t  .  Note  that  the  characteristic 

radius  for  energy  deposition,  R  ,  increases  with  time. 

c 

The  chemical  model  for  these  first  tests  was  taken  to  ba  a  mixture 
of  H2:02:N2  '*"n  the  rat*os  2:1:10  at  1  atm  and  T^=300°K.  The  curve  of  the 
induction  time  as  a  function  of  temperature  for  this  mixture  is  shown  in 
Fig.  (5)  and  was  derived  from  the  detailed  studies  of  the  }i2~°2  mechanisms 
(Burks  and  Oran,  1980) . 

In  the  similarity  model,  k  was  estimated  by  comparing  the  formula 
at  300°K, 


Xs  =  3-2  xlO3  /  T(°K)  erg 

M  52  1  M  Cm  seC  °K'  (27) 

which  assumes  that  an  average  molecular  distance  o  and  an  average 
molecular  weight  M  may  be  found,  to  the  more  exact  formulation 
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where  W.,  is  a  function  of  {X.}  and  the  atomic  masses  (m.)  suggested 
3k  3  3 

by  Mason  and  Saxena.  This  gives  us 


a  =  3.16  A 

M  =  24.3. 

This  approximation  is  valid  because  the  similarity  solution  is  only 
accurate  for  ignition  before  any  major  amount  of  product  or  intermediates 
are  formed.  Then  the  parameter  <  is  found  using  the  definition  in  Eq  (4) . 
Determining  the  temperature  (radius)  at  which  <  is  to  be  evaluated  is  part 
of  the  calibration  to  be  done. 

In  the  first  case  studies, 

R  =  .1  cm 

o 

-4 

T  =  1.0  x  10  sec. 
o 

By  varying  Eq,  the  similarity  solution  indicates  that  the  minimum  ignition 

4 

energy  is  about  3.3  x  10  ergs.  Figs.  (6a  and  b)  show  the  typical  behavior 

4 

of  A(t) ,  R  (t) ,  T (R=0)  and  I(t)  for  3.3  x  10  ergs.  The  functions  I(t)  and 
c 

< (t)  were  evaluated  at  the  central  temperature,  R=0.0cm,  although 
calculations  in  which  these  quantities  were  evaluated  at  R=Rc  were 
virtually  indistinguishable.  Results  from  the  corresponding  detailed 
simulations  are  shown  in  Fig  (7)  for  three  values  of  Eq.  A  comparison 
then  tells  us  that  agreement  between  the  similarity  solution  and  the 
detailed  simulation  is  good.  We  note  here  that  for  this  case  both  models 
predict  ignition  at  essentially  the  same  time  for  a  range  cf  input  energies. 


In  the  second  example. 


R  =  0.025  cm 
o 

-4 

t  =  1  x  10  sec. 
o 

2 

The  similarity  model  predicts  a  minimum  ignition  energy  of  —  8  x  10  ergs. 

4 

However,  for  the  range  of  energies  tested,  up  to  1  x  10  ergs,  the  full 

simulation  does  not  show  ignition,  but  predicts  that  some  burning  does 

occur  and  the  flame  is  eventually  quenched.  Note  however  that  although 

flame  propagation  is  not  predicted  in  these  simulations,  the  values  of 

I(t)  and  <(t)  were  evaluated  at  the  central  temperature.  Evaluating  them 

at  a  further  radius,  R=R^ ,  did  not  improve  agreement.  The  similarity  model 

must  be  calibrated  by  evaluating  both  x(t)  and  I(t)  at  a  particular 

characteristic  radius  which  in  this  case  must  be  even  larger  than  R=R  . 

c 

Conclusion 

We  have  thus  shown  the  existence  of  a  quench  volume  which  has  been 

discussed  previously  in  the  literature  by,  for  example,  Lewis  and  von 

Elbe  (1961) ,  Dixon-Lewis  (1978)  and  Overley  et  at  (1978) .  It  occurs  in 

the  simulation  mentioned  above  because  of  the  relatively  swift  effects 

of  molecular  diffusion  in  the  small  volume  in  which  the  energy  is  deposited. 

We  know  that  for  a  par^iculAr  geometry  and  gas  mixture,  the  quench  volume  must 

be  related  to  R  ,  x  ,  E  and  the  properties  of  the  material  X  (T) ,  {h.(t)}, 
o  o  o  mg 

Y,  and  t  .  We  are  currently  using  the  combination  of  the  similarity  and 
detailed  simulations  to  test  various  approximations  for  the  quenching  volume 
and  to  test  the  sensitivity  of  this  volume  to  the  manner  in  which  chemical 
reaction  is  initiated. 
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Reaction  Rate  Constants 
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Fig.  2  -  Temperature  profile  calculated  from  the  detailed  simulation 
model  for  a  flame  in  an  H2 : O2 :N2/2 : 1 : 10  mixture  ignited  by  a  quadratic 
temperature  profile  at  the  onset  of  the  calculation.  Also  shown  is  a 
graph  of  the  cell  size,  x,  as  a  function  of  position. 


21 


NUMBER  DENSITY  (cm" 


Fig.  3  -  Calculated  profiles  of  species  concentrations  corresponding 
to  the  temperature  profile  in  Fig.  2. 
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Fig.  5  -  Temperature  as  a  function  of  induction  time  for  a  mixture  of 
H2 :(>2 :N2/2 :1 :10  calculated  using  the  reaction  scheme  in  Table  I. 
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